※ 본 문서는 PC 환경에 최적화되어있습니다.

1 서울시 고령자 인구수

1.1 사용한 패키지

library(tidyverse)
library(readr)
library(tmap)
library(tmaptools)
library(maptools)
library(rgdal)
library(DT)
library(spdep)
library(ggformula)
library(sp)
library(spatialreg)
library("classInt")
library(AICcmodavg)

1.2 서울시 고령자 인구수 단계구분도

sigungu_shp <- readOGR(dsn = "./data/shp/Seoul_Dong.shp", encoding = "UTF-8", stringsAsFactors = FALSE)

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Maltese-10\data\shp\Seoul_Dong.shp", layer: "Seoul_Dong"
with 424 features
It has 7 fields

summary(sigungu_shp)

Object of class SpatialPolygonsDataFrame
Coordinates:
       min      max
x 179354.7 216408.6
y 436257.1 466546.9
Is projected: TRUE 
proj4string :
[+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000
+ellps=GRS80 +units=m +no_defs]
Data attributes:
  ADM_DR_CD             동명             Shape_Leng      Shape_Area      
 Length:424         Length:424         Min.   : 2323   Min.   :  215511  
 Class :character   Class :character   1st Qu.: 3931   1st Qu.:  659043  
 Mode  :character   Mode  :character   Median : 4868   Median :  983122  
                                       Mean   : 5495   Mean   : 1431588  
                                       3rd Qu.: 6164   3rd Qu.: 1493250  
                                       Max.   :20854   Max.   :13191610  
   population         objectid       시군구명        
 Min.   :  301.0   Min.   :  1.0   Length:424        
 1st Qu.:  474.2   1st Qu.:106.8   Class :character  
 Median :  612.0   Median :212.5   Mode  :character  
 Mean   : 1536.0   Mean   :212.5                     
 3rd Qu.:  796.0   3rd Qu.:318.2                     
 Max.   :45455.0   Max.   :424.0

old <- read.csv("./data/고령자.csv")
old <- filter(old, !(동 == "소계"))
old <- old[, c(2, 3, 4)]
old2 <- as.character(old$계)
old2 <- gsub(",", "", old2)
old2 <- as.numeric(old2)
old <- cbind(old, old2)
old <- old[, c(1, 2, 4)]
old <- dplyr::rename(old, 동명 = "동")
old3 <- old %>% unite(동명, 자치구, 동명, sep = " ")
map_and_data <- merge(sigungu_shp, old3, by = "동명")
tm_shape(map_and_data) + tm_polygons("old2", id = "동명", convert2density = T, palette = "Blues", 
    border.col = "black", border.alpha = 0.25)

1.3 Moran’s I statistic of spatial autocorrelation

전역적 Moran 1 지수는 특정 공간 단위의 속성 값과 그것의 공간시차와의 상관정도를 나타내 주는 통계치로서 유사한 값들의 전반적인 군집경향을 나타내는 지표이다. 서울시 동별 고령자 수의 공간자기상관을 분석하기 위해 Moran I지수를 살펴보았다.

nb <- poly2nb(map_and_data, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
moran.test(map_and_data$old2, lw, zero.policy = TRUE, na.action = na.omit)


    Moran I test under randomisation

data:  map_and_data$old2  
weights: lw    

Moran I statistic standard deviate = 10.698, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.3044727197     -0.0023640662      0.0008226718

moran.plot(map_and_data$old2, lw, zero.policy = TRUE)

공간가중행렬은 Queen contiguity 방식을 사용하였으며, Moran I 지수는 약 0.3044 으로 산출되었다. 또한 모란 산포도를 그려본 결과 양의 관계가 나타나고 있음을 알 수 있다. 따라서 서울시 동별 고령자 인구분포는 공간적 자기상관성을 가지고 있다고 말할 수 있다. 그러나 이 지역안에서도 발생할 수 있는 공간자기상관의 국지적 변이를 고려하여 LISA 중 국지 Moran을 사용하였다.

1.4 local Moran’s I

로컬 모란 지수를 산출한 결과와 그것을 급간으로 나눠 시각화한 결과이다.

local_m <- localmoran(map_and_data$old2, lw)
datatable(local_m, filter = "top")

local_m <- local_m[, 1]
sigma <- sd(local_m)
mbreak <- c(-Inf, -sigma, 0, sigma, Inf)
group <- cut(local_m, mbreak)
pal <- cm.colors(4)
my.col <- pal[group]
plot(map_and_data, axes = T, col = my.col)

1.5 OLS

basic <- read.csv("./data/수급자.csv")
basic <- filter(basic, !(동 == "소계") & !(동 == "기타"))
basic <- dplyr::rename(basic, 동명 = "동")
basic2 <- as.character(basic$가구)
basic2 <- gsub(",", "", basic2)
basic2 <- as.numeric(basic2)
basic <- cbind(basic, basic2)
basic <- basic[, c(2, 3, 15)]
basic10 <- basic %>% unite(동명, 자치구, 동명, sep = " ")
sp <- merge(basic, old, by = "동명")
gf_point(sp$old2 ~ sp$basic2, xlab = "서울시 기초생활수급자", ylab = "서울시 고령자수") %>% 
    gf_lm(interval = "predict") %>% gf_lm(interval = "confidence")

lm(sp$old2 ~ sp$basic2)


Call:
lm(formula = sp$old2 ~ sp$basic2)

Coefficients:
(Intercept)    sp$basic2  
  19607.283        8.256

summary(lm(sp$old2 ~ sp$basic2))


Call:
lm(formula = sp$old2 ~ sp$basic2)

Residuals:
   Min     1Q Median     3Q    Max 
-21582  -5807   -531   4982  33727 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 19607.283    662.640  29.590  < 2e-16 ***
sp$basic2       8.256      1.106   7.465 4.73e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8779 on 427 degrees of freedom
Multiple R-squared:  0.1154,    Adjusted R-squared:  0.1134 
F-statistic: 55.72 on 1 and 427 DF,  p-value: 4.727e-13

위에서 시행한 회귀분석의 결과 나온 잔차를 시각화하고 잔차의 공간자기상관을 확인하였다.

bfit <- predict(lm(sp$old2 ~ sp$basic2), data.frame(sp$old2))
resi <- (sp$old2 - bfit)
old3 <- cbind(old3, resi)
old_map <- merge(map_and_data, old3, by = "동명")
tm_shape(old_map) + tm_polygons("resi", id = "동명", palette = "Reds", border.col = "black", 
    border.alpha = 0.25)

moran.test(old_map$resi, lw, zero.policy = TRUE, na.action = na.omit)


    Moran I test under randomisation

data:  old_map$resi  
weights: lw    

Moran I statistic standard deviate = 1.863, p-value = 0.03123
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.0510737063     -0.0023640662      0.0008227506

moran.plot(old_map$resi, lw, zero.policy = TRUE)

아주 약한 정도의 자기상관성만을 확인할 수 있다.

1.6 SAR

SAR(spatial autoregressive model)은 종속변수가 공간자기 상관을 갖는 경우로, 공간자기상관의 정보를 독립변수로 회귀식에 설정한다. 현재 회귀모델의 종속변수인 서울시 동별 고령자수의 공간자기 상관은 위에서 확인하였다.

old_map <- merge(old_map, basic10, by = "동명")
mod.sar <- lagsarlm(old_map$old2.x ~ old_map$basic2, data = old_map, listw = lw, 
    zero.policy = T)
summary(mod.sar)


Call:lagsarlm(formula = old_map$old2.x ~ old_map$basic2, data = old_map, 
    listw = lw, zero.policy = T)

Residuals:
     Min       1Q   Median       3Q      Max 
-20545.7  -5229.1  -1129.1   4417.6  32084.3 

Type: lag 
Coefficients: (numerical Hessian approximate standard errors) 
                Estimate Std. Error z value  Pr(>|z|)
(Intercept)    9334.7590  1339.4500  6.9691 3.190e-12
old_map$basic2    6.6077     1.0180  6.4909 8.534e-11

Rho: 0.45502, LR test value: 61.85, p-value: 3.6637e-15
Approximate (numerical Hessian) standard error: 0.052977
    z-value: 8.5891, p-value: < 2.22e-16
Wald statistic: 73.772, p-value: < 2.22e-16

Log likelihood: -4408.147 for lag model
ML residual variance (sigma squared): 63295000, (sigma: 7955.8)
Number of observations: 423 
Number of parameters estimated: 4 
AIC: 8824.3, (AIC for lm: 8884.1)

res <- mod.sar$residuals
sigma <- sd(res)
classes_fx <- classIntervals(res, n = 5, fixedBreaks = c(-Inf, -sigma, 0, sigma, 
    Inf), rtimes = 1)

cols <- findColours(classes_fx, pal)

par(mar = rep(0, 4))

plot(map_and_data, col = cols, border = "black")

legend(x = "bottom", legend = names(attr(cols, "table")), cex = 1, fill = attr(cols, 
    "palette"), bty = "n", ncol = 5, title = "Residuals from SAR Model")

1.7 SEM

SEM(spatial error model)은 오류항의 공간자기상관을 통제하기 위해 활용된다. 오류항의 공간자기 상관성은 위에 제시한 바와 같이 매우 낮은 정도에 해당하였으나, SEM을 시행해보았다.

mod.sem <- errorsarlm(old_map$old2.x ~ old_map$basic2, listw = lw, zero.policy = T)
summary(mod.sem)


Call:errorsarlm(formula = old_map$old2.x ~ old_map$basic2, listw = lw, 
    zero.policy = T)

Residuals:
     Min       1Q   Median       3Q      Max 
-21162.4  -5461.5  -1001.6   4550.8  33092.4 

Type: error 
Coefficients: (asymptotic standard errors) 
                 Estimate Std. Error z value  Pr(>|z|)
(Intercept)    19397.3185   891.9002 21.7483 < 2.2e-16
old_map$basic2     7.2185     1.1143  6.4779 9.302e-11

Lambda: 0.47357, LR test value: 60.857, p-value: 6.1062e-15
Approximate (numerical Hessian) standard error: 0.05501
    z-value: 8.6089, p-value: < 2.22e-16
Wald statistic: 74.113, p-value: < 2.22e-16

Log likelihood: -4408.644 for error model
ML residual variance (sigma squared): 63202000, (sigma: 7950)
Number of observations: 423 
Number of parameters estimated: 4 
AIC: 8825.3, (AIC for lm: 8884.1)

res <- mod.sem$residuals
sigma <- sd(res)
classes_fx <- classIntervals(res, n = 5, fixedBreaks = c(-Inf, -sigma, 0, sigma, 
    Inf), rtimes = 1)
cols <- findColours(classes_fx, pal)
par(mar = rep(0, 4))
plot(map_and_data, col = cols, border = "black")
legend(x = "bottom", legend = names(attr(cols, "table")), cex = 1, fill = attr(cols, 
    "palette"), bty = "n", ncol = 5, title = "Residuals from SEM Model")

1.8 AIC

OLS 모형과, SAR, SEM 모형의 적합도를 AIC(Akaike Information Criterion)를 통해 비교하였다.

ols <- lm(sp$old2 ~ sp$basic2)
AIC(ols)

[1] 9012.198

AIC(mod.sar)

[1] 8824.294

AIC(mod.sem)

[1] 8825.287

비교결과, ols 모형보다 SAR 모형과, SEM 모형의 AIC 값이 더 낮은 것을 확인할수 있다.

2 디트로이트의 911 신고전화

주제를 정하던 도중 신고전화 공간분포의 분석을 해보면 좋을 것 같아 해당 주제를 가지고 분석을 해 보았다.

2.1 사용한 데이터셋

사용한 데이터셋으로는 2016년 9월부터 2018년 1월까지 디트로이트의 911 긴급전화 신고 데이터셋과, 디트로이트가 속한 웨인 카운티와 다른 카운티가 같이 섞여있는 센서스 셰이프파일, 그리고 셰이프의 클립연산으로 디트로이트 부분만 추출하기 위한 디트로이트 영역 셰이프 파일이 있다.

모든 자료는 Data Driven Detroit’s Open Data Portal 및 City of Detroit Open Data Portal 에서 구득하였다.

2.2 사용한 패키지

library(tidyverse)  # 이제는 기본
library(raster)  # 공간정보 처리
library(sp)  #      ''
library(GISTools)  #      ''
library(rgdal)  #      ''
library(tmap)  # 지도 도식
# library(RColorBrewer)
library(spdep)  # 공간통계

2.3 데이터 전처리

최대한 공간적 자기상관이 드러나게끔 신고전화 수를 지역 인구 수로 나눈 것으로 분석을 진행하였다.

2.3.1 미국 센서스 데이터

디트로이트가 속한 Wayne 카운티 및 Oakland, Macomb 카운티가 혼재된 인구 셰이프를 불러왔다. 또한, 셰이프 객체 하나하나에 얼마나 많은 신고 위치 수가 있는지, 클립연산을 행하기 위해 좌표계 통일이 필요하였다. 따라서 해당 좌표계로 통일하기 위해 해당 좌표계를 추출하였다.

spdata <- readOGR("./data/shp", layer = "KEEDPublishingDemographics", encoding = "UTF-8")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Maltese-10\data\shp", layer: "KEEDPublishingDemographics"
with 174 features
It has 15 fields

nad <- crs(spdata)
nad

CRS arguments:
 +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0

2.3.2 911 신고전화 데이터

Municipality(지방자치단체)당 신고건수가 얼마나 있는지 알아보기 위하여 인구 폴리곤과 911 신고위치 데이터를 중첩하여, 폴리곤 객체당 얼마나 많은 수가 있는지 나타내었다.

raw <- read_csv("./data/911_Calls_For_Service.csv")
# raw$call_timestamp <- as.Date(raw$call_timestamp) em2019 <- subset(raw,
# as.Date(raw$call_timestamp) >= as.Date('2018-01-01')) %>% subset(call_timestamp
# <= as.Date('2019-12-31'))

# gcrs <- crs gcrs@projargs <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
# #wgs1984

ecpoints <- SpatialPointsDataFrame(raw[, 24:25], raw, proj4string = nad)


# pcsecpts <- spTransform(ecpoints,crs)

spdata@data$counts911 <- poly.counts(ecpoints, spdata)

tmap_mode("view")
tm_shape(spdata) + tm_polygons("counts911", id = "counts911", alpha = 0.6)

# plot(spdata) plot(ecpoints, add = TRUE)

2.3.3 디트로이트로 지역 국한

폴리곤의 Clip 연산을 수행하여 디트로이트만 나오도록하였다. 구한 디트로이트 시 경계 셰이프는 센서스 셰이프와 CRS를 통일시켜주었다.

bd <- readOGR("./data/shp", layer = "f00e39fd-6902-4d8d-92f7-9f2f5d8156072020329-1-5i4vbz.rjbjg", 
    encoding = "UTF-8")

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Maltese-10\data\shp", layer: "f00e39fd-6902-4d8d-92f7-9f2f5d8156072020329-1-5i4vbz.rjbjg"
with 1 features
It has 11 fields

ndbd <- spTransform(bd, nad)


clip <- raster::intersect(spdata, ndbd)

tmap_mode("view")
tm_shape(clip) + tm_polygons("counts911", id = "NHOOD", alpha = 0.6)

2.3.4 911 신고전화 - 인구 비

인구당 신고전화 수의 비를 알기 위해 신고 수를 인구 수로 나누었다.

c2p_ratio <- vector(mode = "numeric", length(clip@data$GlobalID))
for (i in 1:length(clip@data$GlobalID)) {
    c2p_ratio[i] <- clip@data$counts911[i]/clip@data$Total_Popu[i]
}

c2p_ratio[!is.finite(c2p_ratio)] <- 0
clip@data$c2p_ratio <- c2p_ratio



tm_shape(clip) + tm_polygons("c2p_ratio", id = "NHOOD", alpha = 0.6)

2.4 공간적 자기상관

2.4.1 공간가중행렬

공간적 자기상관을 파악하기 위해 필수적인 공간가중행렬을 구축하였다. 구축방식은 Queen’s Case를 사용하였다.

p2n <- poly2nb(clip, queen = TRUE)  # 인접행렬


plot(clip, col = "gray", border = "blue", lwd = 2)
plot(p2n, coordinates(clip), col = "red", lwd = 2, add = TRUE)

p2n은 폴리곤 데이터 안 인접하는 개체(object)의 목록이다. 예를 들어 첫번째 개체의 인접 목록을 보고 싶다면 다음과 같이…

p2n[[1]]

[1]  9 34

1번 개체의 인접하는 개체 수는 2개가 있다.

clip$OBJECTID.1[1]

[1] 8

clip$OBJECTID.1[c(9, 34)]

[1]  26 129

1번 개체의 부여된 개체코드는 8번이고, 인접하는 두 개체의 코드는 9번과 34번이다.

다음은 인접행렬에 가중치를 부여하여 가중행렬을 만드는 과정이다. 이 사례에서는 인접하는 각각의 개체들에 같은 가중치(1/n)을 부여하였다. 가장 직관적이게 인접 수치를 볼 수 있으나, 연구지역 외곽 개체의 lagged 값에 영향을 주어 잠재적으로 해당 데이터의 진정한 공간자기상관을 overshoot 혹은 undershoot할 수 있다는 단점이 있다.

lw <- nb2listw(p2n, style = "W", zero.policy = TRUE)  # 가중행렬

lw$weights[2]

[[1]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

2.4.2 Spatially lagged values

인접하는 각각의 개체에 가중치를 할당하였다. 즉, 인접긴급전화비율의 평균을 구할 때, 각각의 개체에 해당 가중치가 적용된다.

call.lag <- lag.listw(lw, clip$c2p_ratio)

2.4.3 Moran’s I (Global)

모란 지수 구하기

복잡한 방법

# moran

# 회귀 생성
M <- lm(call.lag ~ clip$c2p_ratio)

Lagged된 비율과 원래 비율을 선형회귀에 피팅하여 나온 회귀선의 기울기가 모란 I 지수이므로 위와같이 회귀 모델을 생성하였다.

# plot
plot(x = clip$c2p_ratio, y = call.lag, main = " Moran Scatterplot Call/Pop Ratio")
abline(lm(call.lag ~ clip$c2p_ratio), lty = 3, lwd = 4, col = "red")

# moran's I index
coef(M)[2]

clip$c2p_ratio 
     0.2588687

산포도와 모란 I 지수이다.

쉬운 방법

혹은 아래와 같은 함수에 집어넣어서 쉽게 구할 수도 있다.

moran.test(clip$c2p_ratio, lw)


    Moran I test under randomisation

data:  clip$c2p_ratio  
weights: lw    

Moran I statistic standard deviate = 2.8384, p-value = 0.002267
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.258868658      -0.020833333       0.009710448

# 함수 plot
moran.plot(clip$c2p_ratio, lw)

유의수준 구하기(Monte Carlo Simulation)

복잡한 방법

몬테 카를로 검정에서는 속성값이 임의로 개체에 지정돼서 모란 값을 계산한다. 출력 결과는 속성 값이 연구 영역 전체에 랜덤하게 분포한다는 귀무 가설에서 모란의 I 값을 샘플링한 분포다. 그런 다음 처음에 구한 모란 I 값을 이 표본 분포와 비교한다. 출력 결과는 다음과 같이 나타내어 p-value를 구할 수 있다.

p-value \(= \dfrac{N_{extreme}+1}{N+1}\)

여기서 \(N_{extreme}\)은 처음에 구한 모란 지수보다 분포의 가장자리에 더 가까운 MC출력 모란지수의 수를 의미하며, \(N\)은 MC 시뮬레이션을 몇 번 돌렸는지를 나타낸다.

###### MC

n <- 599  # 실험 수
I.r <- vector(length = n)


for (i in 1:n) {
    # 신고전화비율을 임의로 섞기
    x <- sample(clip$c2p_ratio, replace = FALSE)
    # lag 값을 새로 구하기
    x.lag <- lag.listw(lw, x)
    # 기울기
    M.r <- lm(x.lag ~ x)
    I.r[i] <- coef(M.r)[2]
}
# Simulated된 모란 지수 도식 우리가 구한 모란 수 또한 빨간 선으로 도식
hist(I.r, main = NULL, xlab = "Moran's I", las = 1)
abline(v = coef(M)[2], col = "red")

# p
N.extreme <- sum(coef(M)[2] < I.r)

p <- (N.extreme + 1)/(n + 1)

p

[1] 0.01

구한 p 값은 지방자치단체 수준에서 군집화 되지 않을 확률이 적다는 점 (p<0.05)을 시사하며, 귀무가설을 기각한다.

쉬운 방법

혹은 아래와 같이 구할 수도 있다.

moran.test(clip$c2p_ratio, lw)


    Moran I test under randomisation

data:  clip$c2p_ratio  
weights: lw    

Moran I statistic standard deviate = 2.8384, p-value = 0.002267
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.258868658      -0.020833333       0.009710448

MC <- moran.mc(clip$c2p_ratio, lw, nsim = 599)

# View results
MC


    Monte-Carlo simulation of Moran I

data:  clip$c2p_ratio 
weights: lw  
number of simulations + 1: 600 

statistic = 0.25887, observed rank = 594, p-value = 0.01
alternative hypothesis: greater

# Plot the distribution
plot(MC, main = "", las = 1)

몬테카를로 시뮬레이션에서 p-value가 다른 이유는 할 때마다 임의로 섞이기 때문이다. 계속된 니트에서도 p 값은 유의수준인 0.05 밑을 유지하였다.

2.4.4 Moran’s I(Local)

국지적 모란지수의 경우 특정 지역의 값과 인접한 주변 지역들이 갖는 값의 가중 평균 값이 서로 유사하게 나타나면 정적인 자기상관으로, 반대로 특정 지역의 값과 인접한 주변 지역들의 가중 평균값과의 차이가 크게 나타나면 부적인 자기상관으로 나타난다.

국지적 모란지수를 아래와 같이 도식하였다.

clip@data$lom <- localmoran(clip$c2p_ratio, lw)[, 1]


tmap_mode("view")
tm_shape(clip) + tm_polygons(col = "lom", title = "Local Moran's I", midpoint = NA, 
    legend.format = list(flag = "+")) + tm_style("col_blind") + tm_scale_bar(width = 0.15) + 
    tm_layout(legend.position = c("left", "bottom"), legend.text.size = 1)

2.4.5 LISA 지표

plot(x = scale(clip$c2p_ratio), y = scale(call.lag), main = " Moran Scatterplot 911calls/pop")
abline(h = 0, v = 0)
abline(lm(call.lag ~ clip$c2p_ratio), lty = 3, lwd = 4, col = "red")

아래 LISA 지표를 만들어보기 위해 정규화를 한 값을 가지고 산포도를 그려보았다.

2.4.5.1 지표 산정

산포도 상에서 점에 해당하는 값을 1 ~ 4사분면 및 유의수준을 기준으로 반복문 및 다중조건문을 사용하여 분류하였다. 여기서 1은 HH, 2는 LL, 3은 HL, 4는 LH, 5는 유의하지 않음이다.

lom <- localmoran(clip$c2p_ratio, lw)

lomp <- lom[, 5]
s.c2p <- scale(clip$c2p_ratio)
s.lag <- scale(call.lag)



# 반복-조건문

flag <- vector(mode = "numeric", length(clip@data$GlobalID))
for (i in 1:length(clip@data$GlobalID)) {
    if (s.c2p[i] >= 0 && s.lag[i] >= 0 && lomp[i] <= 0.05) {
        flag[i] <- 1
    } else if (s.c2p[i] <= 0 && s.lag[i] <= 0 && lomp[i] <= 0.05) {
        flag[i] <- 2
    } else if (s.c2p[i] >= 0 && s.lag[i] <= 0 && lomp[i] <= 0.05) {
        flag[i] <- 3
    } else if (s.c2p[i] <= 0 && s.lag[i] >= 0 && lomp[i] <= 0.05) {
        flag[i] <- 4
    } else {
        flag[i] <- 5
    }
}




clip$flag <- flag

2.4.5.2 도식

여기서는 LISA지표를 통해 국지적으로 공간적 연관성을 분석하였다. 인구대비 신고건 수가 많은 다운타운을 중심으로 High-High 클러스터가 나온 것을 볼 수 있다.

breaks <- seq(1, 5, 1)

labels <- c("high-High", "low-Low", "High-Low", "Low-High", "Not Signif.")

np <- findInterval(clip$flag, breaks)

colors <- c("red", "blue", "lightpink", "skyblue2", "white")
plot(clip, col = colors[np])
mtext("Local Moran's I", cex = 1.5, side = 3, line = 1)
legend("topleft", legend = labels, fill = colors, bty = "n")

2.4.5.3 국지적 모란 I p-value

다운타운 주변을 제외한 나머지 지역이 LISA 지표에서 유의하지 않은 이유가 다음의 국지모란지수의 p 값에서 나타난다.

clip$lomp <- lomp

tm_shape(clip) + tm_polygons(col = "lomp", title = "p-value", breaks = c(0, 0.01, 
    0.05, 0.1, 1), border.col = "black", palette = "-Greens") + tm_scale_bar(width = 0.15) + 
    tm_layout(legend.position = c("left", "bottom"))

2.5 마치며

신고 건수가 많은 다운타운을 중심으로 클러스터가 형성되어 있으며, 여러 지도 도식을 통하여 핫스팟을 파악할 수 있었다.
하지만, 지역을 디트로이트로 한정하면서 수행한 클립 연산시 중첩한 두개의 폴리곤 사이에서 생긴 Sliver를 처리를 하지 못하여, 가중행렬 생성시 폴리곤 가장자리에 있는 작은 개체들까지 가중치 부여가 되었다.
또한 비율을 구할 때, 몇몇 개체의 인구수 속성 값이 0이라 비율을 구하지 못하고 해당 개체는 비율 값을 0이라 두어서 기울기(모란지수)를 구할 때, 영향을 미쳤다는 한계가 존재한다.

3 서울시 자치구의 Cluster 분석

3.1 사용한 패키지

library(ggplot2)
library(rgdal)
library(sp)
library(magrittr)
library(maps)
library(spdep)
library(maptools)
library(ape)
library(spData)
library(classInt)
library(cluster)

3.2 자료입력

sp <- readShapePoly("./data/shp/TL_SCCO_SIG_W")
census <- read.csv("./data/Gustats.csv", header = T)
geocensus <- merge(sp, census, by = "SIG_KOR_NM")

xa = coordinates(sp)[, 1]
ya = coordinates(sp)[, 2]
dis = as.matrix(dist(cbind(xa, ya)))
disi = 1/dis
diag(disi) = 0
r1 = census$alien/census$pop  # 외국인비율
r2 = census$pop/census$hh  # 세대당인구
r3 = census$old/census$pop  # 고령인구비율

3.3 비계층 군집분석(pam) 결과

서울의 25개 자치구 중에서 외국인비율과 고령인구비율을 토대로 5개 그룹으로 군집분석한 결과는 다음과 같다.

geocensus$r1 = geocensus$alien/geocensus$pop
geocensus$r2 = geocensus$pop/geocensus$hh
geocensus$r3 = geocensus$old/geocensus$pop

pc = geocensus[, c("SIG_KOR_NM", "r1", "r3")]
qq = pam(cbind(geocensus$r1, geocensus$r3), k = 5)
pc$rank1 = order(pc$r1)
pc$rank3 = order(pc$r3)
pc$cluster = qq$clustering
plot(pc$r3, pc$r1, type = "n", xlab = "고령인구비율", ylab = "외국인비율", main = "서울 자치구 인구속성")
text(pc$r3, pc$r1, labels = pc$SIG_KOR_NM, cex = 0.6, pos = 3)
points(pc$r3[pc$cluster == "1"], pc$r1[pc$cluster == "1"], pch = 15, col = "red")
points(pc$r3[pc$cluster == "2"], pc$r1[pc$cluster == "2"], pch = 15, col = "blue")
points(pc$r3[pc$cluster == "3"], pc$r1[pc$cluster == "3"], pch = 15, col = "green")
points(pc$r3[pc$cluster == "4"], pc$r1[pc$cluster == "4"], pch = 15, col = "grey")
points(pc$r3[pc$cluster == "5"], pc$r1[pc$cluster == "5"], pch = 15, col = "violet")
legend("topright", c("Cluster no.1", "Cluster no.2", "Cluster no.3", "Cluster no.4", 
    "Cluster no.5"), lty = 1, lwd = 2, col = c("red", "blue", "green", "grey", "violet"), 
    text.col = c("red", "blue", "green", "grey", "violet"))

제1그룹 : 강북구, 도봉구, 은평구, 중랑구, 노원구
제2그룹 : 중구, 종로구, 동대문구
제3그룹 : 광진구, 마포구, 성동구, 관악구, 동작구, 성북구
제4그룹 : 영등포구, 금천구, 구로구, 용산구
제5그룹 : 송파구, 강남구, 서초구, 양천구, 강서구, 강동구

4 서울시 자치구의 공간적 자기상관 분석

4.1 x, y 좌표에 의한 Moran’s I 산정

Moran.I(r1, disi, na.rm = T)

$observed
[1] 0.03064059

$expected
[1] -0.04166667

$sd
[1] 0.03281405

$p.value
[1] 0.02755632

Moran.I(r2, disi, na.rm = T)

$observed
[1] -0.009935369

$expected
[1] -0.04166667

$sd
[1] 0.03138569

$p.value
[1] 0.3120109

Moran.I(r3, disi, na.rm = T)

$observed
[1] 0.1088674

$expected
[1] -0.04166667

$sd
[1] 0.03227261

$p.value
[1] 3.094388e-06

RESULT : 서울 자치구의 x, y 좌표에 의하여 산정한 거리를 토대로, Moran’s I 산정결과는 다음과 같다. 외국인비율과 고령인구비율은 공간적 자기상관이 존재한다. 세대당인구는 공간적 자기상관이 없다.

4.2 삼각망 자료에 의한 Moran’s I 산정

4.2.1 Delaunary 삼각망에 의한 공간적 자기상관

coords = cbind(xa, ya)
ind <- sapply(slot(sp, "polygons"), function(x) slot(x, "ID"))
trinb <- tri2nb(coords)
plot(sp, border = "white", col = "green")
plot(trinb, coords, add = TRUE)
title(main = "Delaunay Triangulation")
text(xa, ya, pos = 3, labels = census$SIG_KOR_NM, cex = 0.6, col = "red")

moran.test(r1, nb2listw(trinb, style = "W"))


    Moran I test under randomisation

data:  r1  
weights: nb2listw(trinb, style = "W")    

Moran I statistic standard deviate = 2.6625, p-value = 0.003879
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.26640488       -0.04166667        0.01338862

moran.test(r2, nb2listw(trinb, style = "W"))


    Moran I test under randomisation

data:  r2  
weights: nb2listw(trinb, style = "W")    

Moran I statistic standard deviate = 0.2804, p-value = 0.3896
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      -0.01067300       -0.04166667        0.01221797

moran.test(r3, nb2listw(trinb, style = "W"))


    Moran I test under randomisation

data:  r3  
weights: nb2listw(trinb, style = "W")    

Moran I statistic standard deviate = 4.6466, p-value = 1.687e-06
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.48688034       -0.04166667        0.01293873

RESULT : 서울 자치구의 Delaunary Triangular Link에 의하여 산정한 Neighbor를 토대로, Moran’s I를 산정한 결과는 다음과 같다. 외국인비율과 고령인구비율은 공간적 자기상관이 존재한다. 세대당인구는 공간적 자기상관이 없다.

4.2.2 KNN(N=4) 삼각망에 의한 공간적 자기상관

knn4 = knearneigh(coords, k = 4)
knn4.nb = knn2nb(knn4)
plot(sp, border = "white", col = "green")
plot(knn4.nb, coords, add = TRUE)
title(main = "KNN : K Nearest Neighbours \n (k=4) Triangular Links")
text(xa, ya, pos = 3, labels = census$SIG_KOR_NM, cex = 0.6, col = "red")

moran.test(r1, nb2listw(knn4.nb, style = "W"))


    Moran I test under randomisation

data:  r1  
weights: nb2listw(knn4.nb, style = "W")    

Moran I statistic standard deviate = 1.9894, p-value = 0.02333
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.19996175       -0.04166667        0.01475247

moran.test(r2, nb2listw(knn4.nb, style = "W"))


    Moran I test under randomisation

data:  r2  
weights: nb2listw(knn4.nb, style = "W")    

Moran I statistic standard deviate = 0.75438, p-value = 0.2253
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.04591695       -0.04166667        0.01347921

moran.test(r3, nb2listw(knn4.nb, style = "W"))


    Moran I test under randomisation

data:  r3  
weights: nb2listw(knn4.nb, style = "W")    

Moran I statistic standard deviate = 4.0471, p-value = 2.592e-05
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.44167741       -0.04166667        0.01426315

RESULT : 서울 자치구의 좌표에서 KNN(k=4)를 적용하여 산정한 Neighbor를 토대로, Moran’s I를 산정한 결과는 다음과 같다. 외국인비율과 고령인구비율은 공간적 자기상관이 존재한다. 세대당인구는 공간적 자기상관이 없다.

이후부터는 고령인구비율을 집중적으로 다루뤄보겠다.

4.3 삼각망 자료에 의한 Local Moran’s I 및 G-통계량 산정

4.3.1 Delaunary 삼각망을 이용한 Local Moran’s I 산정

nb1 <- tri2nb(coords)
lom <- localmoran(r3, nb2listw(nb1, style = "W"))
moran.plot(r3, nb2listw(nb1, style = "W"), pch = 19, xlab = "고령인구비율", labels = as.character(sp$SIG_KOR_NM), 
    main = "Moran Scatter Plot")

moran.plot(as.vector(scale(r3)), nb2listw(nb1), pch = 19, xlab = "고령인구비율", 
    labels = as.character(sp$SIG_KOR_NM), main = "Moran Scatter Plot as Vector")

summary(lom)

       Ii                E.Ii              Var.Ii            Z.Ii        
 Min.   :-0.20772   Min.   :-0.04167   Min.   :0.1007   Min.   :-0.3675  
 1st Qu.: 0.03169   1st Qu.:-0.04167   1st Qu.:0.1237   1st Qu.: 0.1858  
 Median : 0.19441   Median :-0.04167   Median :0.1559   Median : 0.6713  
 Mean   : 0.48688   Mean   :-0.04167   Mean   :0.1653   Mean   : 1.2745  
 3rd Qu.: 0.60626   3rd Qu.:-0.04167   3rd Qu.:0.2042   3rd Qu.: 2.0422  
 Max.   : 2.00084   Max.   :-0.04167   Max.   :0.2847   Max.   : 5.1735  
   Pr(z > 0)        
 Min.   :0.0000001  
 1st Qu.:0.0205659  
 Median :0.2510101  
 Mean   :0.2462355  
 3rd Qu.:0.4263032  
 Max.   :0.6433708

1상한의 강북구는 주변지역과 함꼐 Local Moran’s I 가 높은 값을 나타내는 지역으로 판단된다.

4.3.1.1 Local Moran’s I 분석

pr = data.frame(cbind(lom[, 1], (r3 - mean(r3))/sd(r3)))
plot(pr, xlab = "Moran's I", ylab = "Standardised Rate", xlim = c(-0.5, 2.5), pch = 19)
text(pr[, 1], pr[, 2], pos = 1, sp$SIG_KOR_NM, cex = 0.6)
abline(h = 0, col = "red")
abline(v = 0, col = "red")
title("Plot : Local Moran's I")
mtext(side = 3, "- 서울 자치구 고령인구비율 -")

Local Moran’s I 산정결과는 전체적으로 높은 값을 보이고 있는 것으로 파악된다.

4.3.1.2 Local Moran’s I 그래프

ll = localmoran(r3, nb2listw(nb1, style = "W"))
pa = c("#ffffff", "#f0ff00", "#ffce00", "#ff9a00", "#ff0000")
qk = classIntervals(round(ll, 4), n = 5)
qkp = findColours(qk, pa)
plot(sp, col = qkp)
title("Local Moran's I")
legend("topleft", fill = attr(qkp, "palette"), cex = 0.8, legend = names(attr(qkp, 
    "table")), bty = "n")
text(coords, labels = sp$SIG_KOR_NM, cex = 0.6)
mtext(side = 3, "- 고령인구 비율 -")

4.3.2 KNN 삼각망을 이용한 Local Moran’s I 그래프

비계층적 군집분석 방법인 K-Means(K=5)에 의한 Local Moran’s I에 따른 분포패턴은 다음과 같다.

nb2 <- knn2nb(knn4)
zz = (r3 - mean(r3))/sd(r3)
qq = nb2listw(nb2)
wz = lag.listw(qq, zz)
mlm = lm(wz ~ zz)
aa = mlm$coef[1]
bb = mlm$coef[2]
plot(zz, wz, xlab = "Rate", ylab = "Spatial Lag of Rate")
abline(aa, bb)
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
title(paste("Moran Plot I =", format(round(bb, 4))))
text(zz, wz, pos = 3, sp$SIG_KOR_NM, cex = 0.6)

1상한(주변지역과 같이 높은 값) : 강북구, 도봉구, 중구, 종로구, 은평구, 서대문구, 용산구, 성북구, 노원구
2상한(주변지역에 비해 낮은 값) : 해당 없음
3상한(주변지역과 함께 낮은 값) : 강남구, 송파구, 서초구, 광진구, 양천구, 마포구, 강서구, 가동구, 성동구, 영등포구, 관악구, 금천구, 구로구
4상한(주변지역에 비해 높은 값) : 동대문구, 중랑구, 동작구

4.3.3 Local G 통계량 산정

lgg = localG(r3, nb2listw(nb1, style = "W"))
pa = c("#ffffff", "#f0ff00", "#ffce00", "#ff9a00", "#ff0000")
qk = classIntervals(round(lgg, 4), n = 5, style = "kmeans")
qkp = findColours(qk, pa)
plot(qk, pal = pa, xlab = "고령인구", ylab = "누적확률", main = "누적확률 Plot (K-Means)")

plot(sp, col = qkp)
title("Local G-Statistics : K-Means")
legend("topleft", fill = attr(qkp, "palette"), cex = 0.8, legend = names(attr(qkp, 
    "table")), bty = "n")
text(coords, labels = sp$SIG_KOR_NM, cex = 0.6)
mtext(side = 3, "- 고령인구 비율 -")

한강 이남 vs 한강 이북, 구도심 vs 신도심 간의 명확한 차이가 관찰된다

4.3.4 Local Moran’s I 및 Local G 산정결과 출력

Local Moran’s I 및 Local G 통계량을 출력한 결과는 다음과 같다. 또한 Local Moran’s I의 유의확률이 \(\alpha\)=0.05 이하를 추출하면 도봉구, 은평구, 종로구, 강북구, 강남구, 서초구, 송파구 등 7개 지역이 추출되었다.

aaa = data.frame(round(cbind(lom, lgg), 3), row.names = sp$SIG_KOR_NM)
colnames(aaa) <- c("Moran-I", "E(I)", "V(I)", "Z(I)", "Pr(I)", "G-st.")
print(aaa)

         Moran-I   E(I)  V(I)   Z(I) Pr(I)  G-st.
도봉구     1.838 -0.042 0.285  3.523 0.000  2.303
은평구     0.682 -0.042 0.124  2.057 0.020  2.017
동대문구  -0.118 -0.042 0.156 -0.193 0.577 -0.312
동작구    -0.022 -0.042 0.156  0.049 0.481 -0.874
금천구     0.032 -0.042 0.156  0.186 0.426 -1.012
구로구     0.020 -0.042 0.204  0.138 0.445 -1.213
종로구     1.288 -0.042 0.156  3.368 0.000  2.854
강북구     2.001 -0.042 0.156  5.173 0.000  2.539
중랑구    -0.208 -0.042 0.204 -0.367 0.643 -0.517
강남구     1.558 -0.042 0.204  3.540 0.000 -2.588
강서구     0.184 -0.042 0.204  0.500 0.308 -0.681
중구       0.437 -0.042 0.156  1.211 0.113  1.009
 [ reached 'max' / getOption("max.print") -- omitted 13 rows ]

bbb = aaa[(aaa[, 5] < 0.05), ]
bbb

       Moran-I   E(I)  V(I)  Z(I) Pr(I)  G-st.
도봉구   1.838 -0.042 0.285 3.523 0.000  2.303
은평구   0.682 -0.042 0.124 2.057 0.020  2.017
종로구   1.288 -0.042 0.156 3.368 0.000  2.854
강북구   2.001 -0.042 0.156 5.173 0.000  2.539
강남구   1.558 -0.042 0.204 3.540 0.000 -2.588
서초구   0.606 -0.042 0.101 2.042 0.021 -1.670
송파구   1.590 -0.042 0.204 3.611 0.000 -2.590

공간정보분석

ASS-10 ‘1st Law-abiding Dog’ by: 김하정, 서용훈, 황세민

공간정보분석

1 서울시 고령자 인구수

1.1 사용한 패키지

1.2 서울시 고령자 인구수 단계구분도

1.3 Moran’s I statistic of spatial autocorrelation

1.4 local Moran’s I

1.5 OLS

1.6 SAR

1.7 SEM

1.8 AIC

2 디트로이트의 911 신고전화

2.1 사용한 데이터셋

2.2 사용한 패키지

2.3 데이터 전처리

2.3.1 미국 센서스 데이터

2.3.2 911 신고전화 데이터

2.3.3 디트로이트로 지역 국한

2.3.4 911 신고전화 - 인구 비

2.4 공간적 자기상관

2.4.1 공간가중행렬

2.4.2 Spatially lagged values

2.4.3 Moran’s I (Global)

모란 지수 구하기

복잡한 방법

쉬운 방법

유의수준 구하기(Monte Carlo Simulation)

복잡한 방법

쉬운 방법

2.4.4 Moran’s I(Local)

2.4.5 LISA 지표

2.4.5.1 지표 산정

2.4.5.2 도식

2.4.5.3 국지적 모란 I p-value

2.5 마치며

3 서울시 자치구의 Cluster 분석

3.1 사용한 패키지

3.2 자료입력

3.3 비계층 군집분석(pam) 결과

4 서울시 자치구의 공간적 자기상관 분석

4.1 x, y 좌표에 의한 Moran’s I 산정

4.2 삼각망 자료에 의한 Moran’s I 산정

4.2.1 Delaunary 삼각망에 의한 공간적 자기상관

4.2.2 KNN(N=4) 삼각망에 의한 공간적 자기상관

4.3 삼각망 자료에 의한 Local Moran’s I 및 G-통계량 산정

4.3.1 Delaunary 삼각망을 이용한 Local Moran’s I 산정

4.3.1.1 Local Moran’s I 분석

4.3.1.2 Local Moran’s I 그래프

4.3.2 KNN 삼각망을 이용한 Local Moran’s I 그래프

4.3.3 Local G 통계량 산정

4.3.4 Local Moran’s I 및 Local G 산정결과 출력

ASS-10
‘1st Law-abiding Dog’

by: 김하정, 서용훈, 황세민